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Abstract 

We consider the problem of synthetic aperture radar (SAR) imaging and motion estimation of complex 
scenes. By complex we mean scenes with multiple targets, stationary and in motion. We use the 
usual setup with one moving antenna emitting and receiving signals. We address two challenges: (1) 
the detection of moving targets in the complex scene and (2) the separation of the echoes from the 
stationary targets and those from the moving targets. Such separation allows high resolution imaging 
of the stationary scene and motion estimation with the echoes from the moving targets alone. We show 
that the robust principal component analysis (PC A) method which decomposes a matrix in two parts, 
one low rank and one sparse, can be used for motion detection and data separation. The matrix that 
is decomposed is the pulse and range compressed SAR data indexed by two discrete time variables: the 
slow time, which parametrizes the location of the antenna, and the fast time, which parametrizes the 
echoes received between successive emissions from the antenna. We present an analysis of the rank of the 
data matrix to motivate the use of the robust PC A method. We also show with numerical simulations 
that successful data separation with robust PCA requires proper data windowing. Results of motion 
estimation and imaging with the separated data are presented, as well. 



1 Introduction 

In synthetic aperture radar (SAR), the basic problem [6j [16j [7] is to image the reflectivity supported in a 

x 

set J on the ground surface using measurements obtained with an antenna system mounted on a platform 
flying above it, as illustrated in Figure [I] The antenna emits periodically a probing signal f(t) and records 
the echos D{s,t), indexed by the slow time s of the SAR platform displacement and the fast time t. The 
slow time parametrizes the location r(s) of the platform at the instant it emits the signal, and the fast time 
t parametrizes the echoes received between two consecutive illuminations (0 < t < As). The echoes D(s,t) 
are approximately, and up to a multiplicative factor, a superposition of the emitted signals f(t) time-delayed 
by the round-trip travel-time between the platform r(s) and the locations p of scatterers on the ground 

T(s,p) = 2\v(s)-p\/c. (1.1) 

Here c is the wave speed, assumed constant and equal to the speed of light. 
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Figure 1: Setup for synthetic aperture imaging. 

A SAR image is formed by superposing over a platform trajectory of length (aperture) a the data D(s,t) 
convolved with the time reversed emitted signal (matched filtered), and then backpropagated to points in 

x 

the imaging domain J using travel times. With high bandwidth probing signals and large flight apertures, 
SAR is capable of generating images with roughly ten centimeter resolution at ranges of ten kilometers away 
from the platform. Such resolution cannot be achieved for complex scenes because moving targets appear 
blurred and displaced in the images. Image formation should therefore be done in conjunction with target 
motion estimation. In fact, in applications such as persistent surveillance SAR, tracking and imaging the 
moving targets is one of the primary objectives. 

Because targets may have complicated motion over lengthy data acquisition trajectories, the targets are 
tracked over successive small sub- apertures. Each sub-aperture corresponds to a short time interval over 
which the target is in approximate uniform translational motion. The problem is to estimate this motion 
for each time interval in order to bring the small aperture images of the moving targets into focus. Then, 
the images are superposed to form high resolution images over larger apertures. 

The existing algorithms for motion estimation fall roughly into two categories: The first is for the usual 
SAR setup with a single moving antenna, and the target motion is estimated from the phase modulations of 
the return echoes O O [HI [U EH (24J [13l [15l HH1 [20l [TO]. These algorithms assume that all the targets are in 
the same motion, and are sensitive to the presence of strong stationary targets. The second class of methods 
uses more complex antenna systems [14, 23, 22 , with multiple receiver and/or transmitter antennas. They 
form a collection of images with the echoes measured by each receiver-transmitter pair, and then they use 
the phase variation of the images with respect to the receiver /transmitter offsets, pixel by pixel, to extract 
the target velocity. 

In this paper we consider imaging and motion estimation of complex scenes with the usual SAR setup, 
using a single antenna. We address two challenges: (1) the detection of moving targets in the scene and (2) 
the separation of data in subsets of echoes from the stationary scene and echoes from the moving targets. 
The stationary scene can be imaged by itself after such separation, and the motion estimation can be carried 
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out on the echoes from the moving targets alone. We propose and analyze a detection and data separation 
approach based on the robust principle component analysis (robust PC A) method [4 j. Robust PC A is 
designed to decompose a matrix into a low rank one plus a sparse one. The main contribution of this paper 
is to show with analysis and numerical simulations that by appropriately pre-processing and windowing the 
SAR data we can decompose it into a low rank part, corresponding to the stationary scene, and a sparse 
part, corresponding to the moving targets. Our theoretical and numerical study describes the rank of the 
pre-processed SAR data as a function of the velocity, location, and density of the scatterers. It specifies 
in particular how slowly a target can move and still be distinguishable from the stationary scene. It also 
addresses the question of proper windowing of the data for the separation with robust PC A to work. 

The paper is organized as follows: We begin in section [2] with a brief description of basic SAR data 
processing and image formation. There are two processing steps that are key to the data decomposition: 
pulse compression and range compression. We illustrate with numerical simulations in section [3] that robust 
PCA can be used for motion detection and data separation if it is complemented with proper data windowing. 
The analysis is in section [4] Additional numerical results on data separation, motion estimation and imaging 
are in section [5] We end with a summary in section [6j 

2 Basic SAR data processing and image formation 

The antenna emits signals that consist of a base-band waveform /#(£) modulated by a carrier frequency 
v = ^ /(2tt), 

f(t) = co8(u> t)f B {t). (2.2) 

Its Fourier transform is 

f(uj) = Jdtf(t)e^ = ^[f B (uj^uj )^f B (uj-uj )\, (2.3) 

with /#(cj) supported in the interval [—tvB : tvB], where B is the bandwidth. The support of f(u) is the 
same interval with center shifted at cj , and its mirror image in the negative frequencies. 

Imaging relies on accurate estimation of travel times. Recall that the echoes are approximately, up to 
some amplitude factors, superpositions of the emitted signals delayed by the round trip travel times between 
the antenna and the targets in the scene. Suppose that 

f B (t) = <p(Bt), 

with ip a function of a dimensionless argument and compact support in the unit interval [0,1]. Then /#(£) 
has the form of a pulse, and the travel times can be estimated with precision 1/5, the pulse support. But 
such pulses are almost never used in SAR because of power limitations at the antenna [T6 . The SAR echoes 
should be above the antenna's noise level, so the emitted signals should carry large power. However, the 
instantaneous power at the antenna is limited, while large net power can be delivered by longer signals, such 
as chirps. The problem is that the received scattered energy is spread out over the long time support of 
such signals, making it impossible to resolve the time of arrival of different echoes. This is overcome by 
compressing the long echoes, as if they were created by an incident pulse. This pulse compression amounts 
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Figure 2: Left: Pulse compressed data D p (s,t) = f(t — r(s,p)) from a single scatterer at p = (0,5,0), 
recorded by a SAR antenna moving at speed V = 70 m/s along a linear aperture of 400 m, at range L = 10 
km from p. Note the hyperbolic curve defined by the peak. Right: Pulse and range-compressed data D r (s, t') 
with respect to point p Q = (0, 0, 0). The abscissa is the shifted fast time t r . Note that most of the slow time 
dependence is removed in the range-compressed data. 

to convolving the data with the complex conjugate of the time reversed emitted signal f(—t) 

D p (s,t) = J dt'D(s,t')f(t'-t). (2.4) 
The result is as if the antenna emitted the signal 

/„(*) = J dt'f(t')f(t'-t), 

which turns out to be a pulse of support 1/5, as explained in detail in [T6] . 

Another common data pre-processing step is range compression. It amounts to evaluating the pulse 
compressed data at times offset by the round trip travel time to a reference point p Q G J Z ', 

D r (s,lf) = D p (s,lf + T(s,p )). (2.5) 

Here t' is the fast time shifted by r(s,p ), so that 

t' + r(s,p o )=te[0,As}. 

In the Fourier domain, the range compression amounts to removing the large phase ur(s,p ) of D p (s,uj). 
This is obviously advantageous from the computational point of view. But range compression plays a bigger 
role in our context. It is essential for the robust PC A method to work. To illustrate this, we show in Figure 
[2] the pulse compressed echo from a single scatterer, for a linear aperture. The scatterer is at p = (0,5,0), 
with the first two coordinates defining the location in the imaging plane and the third the elevation, which 
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is always zero. The coordinates in the imaging plane are range and cross-range, with origin at the reference 
point. The range is the coordinate of p along the direction pointing from the SAR platform (at the center 
of the aperture) to p Q . The cross-range is the coordinate of p in the direction orthogonal to the range. 

We plot in Figure [2] the amplitude of D p (s,t) as a function of s and t, and note that the location of the 
peak, defined by equation 

f) 2 = |r»-p| 2 , 

is a hyperbola for the linear aperture. The peak lies on some other curve in the (s, t) plane for other 
apertures. The amplitude of the range compressed data is shown in the right plot of Figure [2] We note that 
the dependence on the slow time s has been approximately removed by the range compression. Explicitly, 
it lies on the curve in the (s, t') plane defined by equation 

ct' 

— = \ f(s)-p\ - \r(s)-p \. 

This curve is close to the vertical axis t' = because p and p Q are close to each other, 

\p-p \ <C \r(s)-p\ : Vs. 

Consequently, the matrix with entries D r (s,t), sampled at discrete s and t, appears to be of low rank and 
it can be handled by the robust PCA method. 

We work with the pulse and range compressed data from now on, and to simplify notation, we drop the 
prime from the shifted fast time t' . We borrow terminology from the geophysics literature and call the pulse 
and range compressed echoes data traces. The image is formed by superposing the traces over the aperture, 
and backpropagating them to the imaging points p G J using travel times 



n/2 

D r (sj,T(sj,p X ) -r(sj,p )) 

j=-n/2 

1 f S(a) X 

— dsD r (s,r(s,p ) -r(s,p )). (2.6) 

As J-S(a) 



Here Sj are the discrete slow time samples in the interval [—5(a), S (a)] defining the aperture a along the 
flight track. The sampling is uniform, at intervals As, and 

2S(a) = nAs, with n even. 



Assuming a large n, that is a small As, we approximate the sum in (2.6) by an integral over the aperture 



When the aperture is very large, the data in (2.6 ) is weighted by a factor that compensates for geometrical 
spreading effects over the long flight track, and thus improves the focus of the image. Here we work with 
small apertures where geometrical spreading plays no role, which is why there are no weights in the imaging 



function (2.6). 
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3 Robust PCA for motion detection and SAR data separation 



We begin in section |3.1| with a brief discussion of the robust PCA method. Then, we give in section |3.2| 
a heuristic explanation of why it makes sense to use it for motion detection and data separation. We also 



illustrate in section 3.3 the difficulties arising in the separation, and the improvements achieved by proper 
data windowing. 



3.1 Robust PCA 

The robust PCA method, introduced and analyzed in [4 , applies to matrices M G R niXU2 that are sums 

of a low rank matrix C Q and a sparse matrix S Q . It solves a convex optimization problem called principle 
component pursuit: 

min ||£||*+T7||S||i (3.1) 

£,5ei ri i Xn 2 

subject to C + S = M, (3.2) 

with 

v= i ) x - (3 ' 3) 

v / max{ni,n 2 } 

Here ||£||* is the nuclear norm, i.e. the sum of the singular values of £, and ||<S||i is the matrix 1-norm of S. 
The optimization can be done for any matrix, but the point is that if M = C Q + with C Q low rank and S Q 
sparse and high rank, then the principle component pursuit recovers exactly C Q and S Q . The analysis in [4 
gives sufficient conditions under which the decomposition is exact. These conditions are bounds on the rank 
of C and the number of non-zero entries in the high rank matrix S Q . They are not necessary conditions, 
meaning that the decomposition can be achieved for a much larger class of matrices than those fulfilling the 
assumptions of the theorems in [4]. This is already pointed out in [4]. 



3.2 The structure of the matrix of range compressed SAR data 

We illustrate in this section the pulse and range compressed echoes (the traces) from a complex scene. The 
point of the illustration is to show that typically, the sampled traces from the stationary scene form a low 
rank matrix, whereas those from moving targets give a high rank but sparse matrix. We refer to section [5] 
for the description of the setup of the numerical simulations used to produce the results presented here and 
in the next section. 

Let M be the matrix with entries given by the data traces sampled at the discrete slow and fast times 



where 



and 



Mji =D r ( Sj ,ti), 

Sj = jAs, j = -n/2,...,n/2, 
U I At, I -ra/2, . . . ,m/2. 



(3.4) 

(3.5) 
(3.6) 
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Figure 3: Desired separation of the range and pulse compressed SAR data. Left: \D r (s,t)\ for a sample 
scene with thirty stationary targets. The vertical traces correspond to the stationary targets, while the 
sloped trace corresponds to the moving target. Middle: The stationary component of \D r (s,t)\. Right: The 
moving component of \D r (s,t)\. We plot the absolute value of these matrices only to increase contrast for 
the figures. 

The sampling is at uniform intervals As and At, and n and m are positive integers satisfying 

2S nAs, and As = mAt. (3.7) 

We assume throughout that n and m are even and large. 

We show in the left plot of Figure [3] the matrix M for a scene with a single moving target and thirty 
stationary scatterers. The ideal separation of this matrix would be 

M = C + S 

where C Q is given by the echoes from the stationary targets alone and S Q is given by the echo from the 
moving target. We plot C Q and S Q in the middle and right plot of Figure [3] It appears to the eye that C Q is 
a matrix with almost parallel columns, so we expect it to be low rank. The matrix S Q is given by the sloped 
curve in the (s, t) plane. The range compression does not remove the s dependence of the echo from the 
moving target, as is the case for the stationary ones, and this is why we see this sloped curve. Consequently, 
S has higher rank than C Q . But it is sparse, because we have only one moving target, i.e., one sloped curve. 
Obviously, the matrix S Q will remain sparse for a small number of moving targets, as well. 

Thus, the matrix M appears to have the appropriate structure for the robust PCA method to be useful. 
But we show next that robust PCA alone will not do the job adequately. It must be complemented with 
proper data windowing. 

3.3 Data windowing for separation with robust PCA 

We illustrate with two numerical simulations that robust PCA cannot be applied as a black box to the SAR 
data traces and produce the desired separation of the stationary and moving target parts. However, if it is 
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Figure 4: Robust PCA separation of traces from a scene with seven stationary targets and a moving target. 
From left to right we show the matrix M of traces, the low rank £, and sparse S components from robust 
PCA. Top row is robust PCA applied to the entire data, bottom row is robust PCA applied to a windowed 
part of the data. In each plot the matrices are normalized by the largest value of \D r {s,t)\. The sparse 
components shown in the right column are plotted in (decibel) dB scale to make the contrast more visible. 



applied to properly calibrated subsets (or windows) of the data, robust PCA can separate the traces for a 
large class of complex scenes. 

The results of the first simulation are in Figure [4j We have a scene with seven stationary targets at 
coordinates (0,0, 0)m, (±5,0,0)m, (0,±5,0)m and (±10,0,0)m in the imaging plane, which is at elevation 
zero. There is also a moving target that is at location (0, 0, 0), at slow time s = corresponding to the center 
of the aperture. The target moves with velocity u = ^=(1, l)m/s in the plane. We refer to section jHJ for a 
detailed description of the setup for the simulations. In the top row of plots in Figure [4] we show the results 
with robust PCA applied to the matrix M of data traces shown on the left. The middle and right plots are 
the resulting low rank and sparse parts given by the robust PCA algorithm. Due to the fact that the entire 
matrix M is sparse to begin with, much of the stationary target data is captured by the sparse component. 
Thus, we do not have a good separation. However, the result is very good when we apply the robust PCA 
algorithm on a smaller fast-time window of the data, as shown in the second row of plots. The number of 
nonzero entries in the windowed M is no longer small relative to its size, and robust PCA separates the 
moving target trace from the rest. 

In Figure [6] we show the results for a scene with thirty stationary targets and a moving one. The data 
traces from this scene are displayed in Figure [5] The moving target is as in the previous example. The 
point of the simulation is to show that when we work in a large fast time window the traces from all the 
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Figure 5: SAR scene with thirty stationary targets and a moving one. 



stationary targets may no longer form a low rank matrix. The top row of plots in Figure [6] shows that the 
sparse component of the matrix M, as returned by the robust PCA algorithm, has a large residual part from 
the stationary targets. The much improved results in the bottom are obtained by applying robust PCA 
on successive small fast time windows of the data, and then reassembling the separated traces. To give an 
idea of the size of the windows, the matrix M has dimensions n — 296 by m = 16, 384. The matrices are 
windowed in fast-time and are of size 296 x 450. 

The two examples given above show that data windowing plays an important role in achieving a successful 
data separation with robust PCA. The last simulation also illustrates the role of robust PCA in the detection 
of the moving target. While the trace from this target is faint and difficult to distinguish from the others 
over most of the fast time interval in Figure [5| it is clearly visible in Figure [6] after the processing of the 
traces with robust PCA. 



4 Analysis of rank of the matrix of SAR data traces 

We present here an analysis of the rank of the matrix of SAR data traces for simple scenes with one or two 
targets. The goal is to understand how the rank depends on the position of the targets and their velocity. 
We limit the analysis to at most two targets to get a simple structure of the matrix MM T , for which we can 
calculate the rank almost explicitly. The numerical results presented above and in section [5] show that the 
data separation works for complex scenes, with many stationary targets. 



We begin in section 4.1 with the scaling regime used in the analysis. We illustrate it with the GOTCHA 



Volumetric SAR data set [5 for X-band surveillance SAR. The model of the matrix M is given in section 



4.2 We analyze its rank in section 4.3 for a single target, and in section 4.4 for two targets. We end with a 



brief discussion in section I3~5 
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Figure 6: Robust PCA separation of traces from a scene with thirty stationary targets and a moving one. 
The top row of plots shows the results with robust PCA applied to the entire data matrix in Figure |5j The 
low rank and sparse components are shown in the middle and on the right. The bottom row of plots shows 
the results with robust PCA applied to successive small time windows of the data. The sparse components 
shown in the right column are plotted in dB scale to make the contrast more visible. 



4.1 Scaling regime and illustration with GOTCHA volumetric SAR 

The important scales in the problem are: The central frequency v Ql the bandwidth B, the typical distance 
(range) L from the SAR platform to the targets, the aperture a, the magnitude |u| of the speed of the 
targets, the speed V of the SAR platform and i? x , the diameter of the imaging set J . We assume that 
they are ordered as follows 

B < i/ OJ (4.1) 



and 



and we let 



a<L, R <L, 



L = \m-PoV 

We also suppose that the target speed is smaller than that of the SAR platform 

lul < V. 



(4.2) 
(4.3) 

(4.4) 



As an illustration, consider the setup in GOTCHA Volumetric SAR. The central frequency of the probing 
signal is vq = 9.6GHz and the bandwidth is B = 622MHz, so assumption (4.1) holds. The SAR platform 



trajectory is circular, at height H — 7.3km, with radius R = 7.1km and speed V = 250km/h or 70m/s. One 
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circular degree of trajectory is 124m. The pulse repetition rate is 117 per degree, which means that a pulse 
is sent every 1.05m, and As = 0.015s. A typical distance to a target is L = 10km and we consider imaging 



domains of radius R of at most 50m, so assumptions (4.2) hold. The target speed is |u| ~ lOOkm/h or 



28m/s, so it satisfies (4.4). 



For a stationary target we obtain from basic resolution theory (for single small aperture imaging) that 
the range can be estimated with precision c/B = 48cm, and the cross range resolution is XqL/cl = 2.5m, 
with one degree aperture a and central wavelength Ao = 3cm. The image of a moving target is out of focus 
unless we estimate its velocity and compensate for the motion in the imaging function. 

4.2 Data model 

Our model of the data assumes that the scatterers lying on the imaging surface behave like point targets. 
We neglect any interaction between the scatterers, meaning that we make the single scattering, Born ap- 
proximation. The pulse and range compressed data is approximated by 

N ( v 

Dr ^ t)K ^ (4nm-t q (sW fp{t ~ (T(S ' U 8)) ~ T(S ' ^ ))} ' (45) 
in the case of N targets at locations 

Pq( S ) = (Pq( S )^)i Q=1,'-,N, 

with reflectivity a q (u ). Here we used the so-called stop- start approximation which neglects the displacement 
of the targets during the round trip travel time. This is justified in radar because the waves travel at the 
speed of light that is many orders of magnitude larger than the speed of the targets. We refer to [2 for a 



derivation of the model (4.5) in the scaling regime described above. 



Because of our scaling assumptions (4.2), we can approximate the amplitude factors as 

1 1 



47r|r(s) - p q (s)\ 47rL' 

and obtain the simpler model 

/ 1 \ 2 N 

Dr{s ' t)= {^Lj E^(*- Ar ^^W)), (4-6) 

^ ^ q=l 

where we let 

Ar(s, p(s)) = r(s, p(s)) - r(s, p ). (4.7) 
We assume henceforth, for simplicity, that the targets are identical 

a q = a, g = l,...,iV, 



11 



and write 



with 



D r {s,t) K £^M(s,t), 



(4ttL) s 



N 



M(s,t) = ^2f p (t-Ar(s,p q (s))). 

9=1 



The matrix of traces analyzed below is given by discrete samples of (4.9), 



Mji = M (sj_»_i,i;_i) , j = l,...,ra+ 1, Z = l,...,m+1, 



(4.8) 



(4.9) 



(4.10) 



with slow times Sj and fast times defined in (3.5) and (3.6). We also take for convenience a compressed 



pulse given by a Gaussian modulated by a cosine at the central frequency, 

f p (t) = cos(uj t)e- BH ^ 2 . 

4.3 Analysis of rank of the data traces for one target 



(4.11) 



In the case of one target at location p(s) = (p(s),0), the entries of matrix (4.10) are given by 

B 2 



Mji = cos [uj (t — Ar(s, p(s)))] exp 



-(t-Ar(s,p(s))y 



(4.12) 



with At defined by (4.7). The target is moving at speed u = (u, 0), so we have that 

p(s) = p + su, \s\ <5(a), (4.13) 

where we let 

P ■= 0(0) 

be the location of the target at the time s = 0, corresponding to the center of the aperture. 

We study the rank of M, which is equivalent to studying the rank of the symmetric, square matrix 
C e R(^+i)x(^+i) ? w ith entries given by 



C hi = C (^-f-i,s/-f-i) , j,l = l,...,n + l, 



(4.14) 



in terms of the function 



m /2 „oo 

C(s, s') = D ^ s ' t q )Dr(s', t q ) « — / dt 2? r (5, t)2? r (5, t). (4.15) 

q— — m/2 ' 00 



Here we used the assumption that At is small enough to approximate the Riemann sum over q by the integral 
over t. Because the traces D r {s,t) vanish for \t\ > As/2, we extended the integral to the whole real line. 

We obtain after a calculation given in appendix [A] that if the aperture a is small enough, the matrix C 
has a Toeplitz structure. 
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Proposition 1 Assuming that the Fresnel number is bounded by 

a 2 . f L V 



A ( ,L <<min |^'M}' (4 ' 16) 



the matrix C is approximately Toeplitz, with entries given by 



■\/lT 

CO, s') « cos [u; a(s - s')] exp 



The dimensionless parameter a depends linearly on the velocity in the range direction and the cross-range 
offset of the target with respect to the reference point p Q . It is given by 

2um Q 2Vt-F (p-p ) 2u-P (p-p ) 

a ~— cL + cL ' {4 ' 18) 

with m Q the unit vector pointing in the range direction from the center r(0) of the aperture 

m °-|?(0)-p o |' (419) 

P = /-m m^. (4.20) 



and F Q the orthogonal projection 
The unit vector t is defined by 



ds 

It is tangential to the flight track, at the center of the aperture 



dm - Vt. (4.21) 




0.5 1 1.5 2 2.5 3 3.5 4 

Slow time (seconds) 



Figure 7: The matrix C for one stationary target at p = (0, 15,0)m. The reference point is p Q = 0. The 
aperture is 310m, the range L = 10km and the central wavelength is A = 3cm. The plot shows that the 
matrix is essentially constant along the diagonals, it is approximately Toeplitz. 

As an illustration, we plot in Figure [7] the matrix C for an aperture a = 310m, in the GOTCHA regime 
with central wavelength A G = 3cm and range L = 10km. We have a stationary target at p = (0, 15, 0)m, and 
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Figure 8: Matrix M of data traces for a stationary target located at p = p Q = (0,0, 0)m (left) and at 
p = (-10, -10,0)m (right). 



p Q is at the origin. The Fresnel number 



XnL 



320.3 



(4.22) 



is only half the ratio L/R , with R = 15m, and yet the matrix C is essentially Toeplitz. 
Remark: Since 

11^)11! = ^, v-, 

the angle between two rows of the matrix M of data traces, indexed by s and s f , is given by 



cos Z(D r (s-), D r (s' , •)) « cos[cj a(s — s')) exp 



(Ba) 2 (s-s f ) 



f\2 



(4.23) 



Nearby rows with indices satisfying 



a\s 



< 



B' 



are nearly parallel, with a small, rapidly fluctuating angle between them. But rows that are far apart, with 
indices satisfying 

3\/2 



a\s 



> 



B 



are essentially orthogonal, because the right hand side in (4.23) is approximately zero. When the target is 
stationary, and its cross-range offset is zero, then a = 0, and all the entries in C are constant and equal to 
one. Moreover, all the rows of M are parallel to each other, and the matrix has rank one. When we view 
the data traces we see a vertical line at t = 0, as in the left plot of Figure [8j Obviously, the larger \a\ is, the 
closer the indices of the rows that are nearly orthogonal. So we expect the rank of M to increase with \a\. 
This is indeed the case, as we show next. We also illustrate this fact in the right plot of Figure [5J where we 
show the matrix of data traces for a stationary target that is offset from the reference point p Q . Note that 



by definition (4.18), \a\ is large when the cross-range offset of the target and/or the target speed are large. 
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4.3.1 Asymptotic characterization of the rank 



Because matrix C is Toeplitz and large, we can use the asymptotic Szego theory [T7J [3] to estimate its rank. 
To do so, let us define the sequence {cj}j e z with entries 



Cj = ^ A . e 4 cos(£j), £ = B\a\As. 



(4.24) 



3 2BAf 

A finite set of this sequence, for indices \j\ < n, defines approximately the diagonals of the matrix C, 

Cjjncj-t, = l,...,n+l. (4.25) 

Since multiplication of a large Toeplitz matrix with a vector is approximately a convolution, and since 
convolutions are diagonalized by the Fourier transform, it is not surprising that the spectrum of C is defined 
in terms of the symbol Q(0), the coefficients of the Fourier series of (4.24), 

oo 

Q{9)= E c i e ^> 0e(-7r,7r). 

j=-oo 



(4.26) 



With this symbol, we can characterize asymptotically in the limit n — >• oo the rank of C, using Szego's first 
limit theorem [3], that gives 



A/"(n;/3i,/3 2 ) = 

2tt 



lim 

n-^oo n + 1 



l WlM (Q(0))dO. 



(4.27) 



Here l[/3 1? /3 2 ] is the indicator function of the interval /fe] and A/"(n; /3i, fc) is the number of eigenvalues of 
C that lie in this interval. 

We show in appendix IB] that the symbol is given approximately by 



Q(0) 



where 7 G (— 7r, 7r) is defined by 



2BAt£ 



, exp 



r (^-7) 2 i 




r (^+t) 2 i 


[ e 2 J 


+ exp 


L c 2 J 



}' 



7 = [(cjq^As + 7r) mod 2tt] — tt. 



We use this result and (4.27) to obtain an asymptotic estimate of the essential rank, defined by 



rank[C] := A* (n; e||Q|| IIQII c 



sup |Q(0)|. 

9G(-7T,7r) 



(4.28) 



(4.29) 



(4.30) 



Here < e <C 1 is a small threshold parameter, and ||Q||oo is of the order of the largest singular value of 
C. It follows from the Szego theory [TTJ |3] that this singular value is given by the maximum of the symbol, 
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Estimated Rank 
Computed Rank 



5 10 15 20 25 30 5 10 15 

Cross-Range Range Velocity 

(a) Stationary (b) Moving 

Figure 9: Comparison of the computed and estimated rank of C for a single target. Left: Stationary target 
at various cross-range positions. Right: Moving target with various velocities. 



which is of the order of 7r/(2BAt£ > ). We obtain that for n ^> 1 

rank [C] 



n + 1 h J\\\Q\\^oo)(Q(0))de 

'2\a\BAsJ\ogl/e \ 

v 6 1 , 1 , (4.31) 



mm 



where the last equality follows from direct calculation. 

As an illustration, we show with green in Figure |9^i the computed rank of the matrix C for a single 
stationary target, versus the cross-range position. Positions with larger cross-range result in larger \a\ and 



thus larger rank, as expected. The asymptotic estimate (4.31) of the rank is shown in blue. It exhibits 
the same linear growth in cross-range and therefore in |a|, but it is lower than the computed rank. This is 
because in this simulation n is not sufficiently large. We get a good approximation when Cj « for j w n, 



so we can approximate the series (4.26) that defines the symbol by the truncated sum for indices \j\ < n. 



This is not the case in this simulation, so there is a discrepancy in the estimated rank. However, the result 



improves when we increase n, by increasing the aperture. This is illustrated in Figure 10 , where we show 



the rank normalized by the size of the matrix, and note that it has the predicted asymptote as n increases. 



The plot in Figure [9p shows the computed rank (in green) and the asymptotic estimate (4.31) (in blue) 
for a moving target located at p = p Q at s = 0. We plot the rank as a function of the range component of the 
velocity. The rank increases with the velocity, and therefore with |a|, as expected. Moreover, the asymptotic 
estimate is very close to the computed one, because in this case the entries in the sequence {cj}j e z decay 
faster with j. Finally, notice that even for small velocities, the rank is much larger for the moving target 
than the stationary target. 
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20 40 60 80 100 120 

Aperture Length (seconds) 

Figure 10: Convergence of the rank of C normalized by the size (n + 1). The blue line is the computed value 
and the green line is the asymptotic estimate. 

4.4 Analysis of the rank of the data traces for two targets 

In the case of two targets at locations pj(s) for j = 1, 2, the entries of the matrix M follow from equations 



(4.9) and (4.10) 



Mji = ^cos [u (t - Ar(s,pj(s)))] exp 
We take for simplicity the case of two stationary targets 



B 2 



(t-Ar{s,pAs))f 



(4.32) 



p J (s) = p J :=p J (0), J = 1,2. 
Extensions to moving targets with speeds Ui and U2 are straightforward. They amount to redefining the 



parameters ol\ and &2 defined below by adding two linear terms in the target velocity, as in equation (4.18). 

The expression of matrix C = MM T is given in the following proposition. It is obtained with a calculation 
that is similar to that in appendix [Aj using the same assumption on the Fresnel number as in Proposition [I] 



Proposition 2 Assume that the Fresnel number satisfies the bound with u = since the targets are 

stationary. The matrix C has entries defined by the function 



C(s,s f ) 



2BAt 



\ J2cos[u aj(s - s f )] exp 



(Baj) 2 (s — s') 2 



■ cos[co> (ai<s — + /?)] exp 
- cos[(j0 o (aiS r — (%2S + j3)) exp 



B 2 (aiS — QL2S' 



B 2 {a 1 s' -ol 2 s + PY 



(4.33) 
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sampled at the discrete slow times. Here we let 



2Vt-¥ (p j -p ) 



cL 



and 



[m Q • (jij - p )} 2 
2L 



(4.34) 



(4.35) 



The structure of the matrix C is now more complicated. The first term in (4.33) gives a Toeplitz matrix, 
as before. The other two, which are due to the interaction between the targets, give matrices that are 
approximately g- Toeplitz or g-Hankel, depending on the sign of the ratio a^/c^i. 

Definition 1 An (n + 1) x (n + 1) g-Hankel matrix H with shift g G Z + is defined by a sequence {hj}j e ^ as 

Hji = hj-i+g(i-i)i j, / = 1, . . . , n + 1. 

The matrix is Hankel when g = 1. A g- Toeplitz matrix is defined similarly, by replacing g with —g. 
To analyze the spectrum of matrix C, we choose the reference point p Q in such a way that 



a 2 

0L\ 



< 0, and g := 



OL2 



GN. 



Then, C is given by 
with Toeplitz matrix 
and g-Hankel matrix 



C = T + H + H 7 



T ji = c 3-u j> 1 = !,••• ,n + 1, 



Hji = h(j-i)+ g (i-\), j, I = 1, . . . ,n + 1. 



Here {cj}j e z and {/^j}jeN are sequences with entries 



e 4 cos(4ij) + e 4 cos(4 2 j) 



and 



hi 



2BAt 



e 4 cosgiO + C)], C 



6 = %|As, £ = 1,2 



lailAs" 



(4.36) 

(4.37) 
(4.38) 
(4.39) 

(4.40) 
(4.41) 



An example of matrix C and its Toeplitz and g-Hankel parts is shown in Figure 11 
If the range offsets are large, so that 



m = B\p\ 



2B 



£(-l) J {m -(p i -p ) + — 



2L 



>1, 



the g-Hankel matrix has small entries, and C is approximately Toeplitz, as in the single target case. The 
difference of the travel times between the SAR platform and such targets is larger than the compressed pulse 
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Toeplitz Component g-Hankel Component g-Hankel Transpose Matrix C 




01234 01234 01234 01234 

Slow time (seconds) Slow time (seconds) Slow time (seconds) Slow time (seconds) 



Figure 11: Components of matrix C for two stationary targets located at p x — (0.15, 15, 0)m and p 2 = 
(—0.15, — 5,0)m. On the left we show the Toeplitz part T. The next two plots show the g-Hankel part H 
and its transpose H T . The right plot shows the sum, i.e., the matrix C. 



width, and their interaction in (4.33) is negligible. If the range offsets are small, the structure of the matrix 



C is as in equation (4.37), and the estimate of its rank follows from the recent results in [I9j [TTJ [12] . They 
say that the g-Hankel terms H + H T have a negligible effect on the rank in the limit n — » oo. See appendix 
[C]for more details. Thus, in either case, the rank estimate of C is given by equation fl4.27b, m terms of the 



symbol Q(&) defined by (4.26), using the sequence {cj}j e z with entries (4.40). 



Explicitly, the symbol is given by 



Q(0) 



2BAt^ 



.fail: 



.(fayiT 



2BAt£ 2 



_ (fa72T 



(4.42) 



with 7j G (— 7r,7r) defined by 



7j = [(uJoajAs + 7r) mod 2tt] — 7r, 



for j = 1, 2, and the rank is given by 



rank [C] 



(n + 1) r 



2tt 



l[e\ moo ,oo)(Q(0))d0. 

J — 7T 



(4.43) 



(4.44) 



4.4.1 Illustration 



We compare in Figure 13 the computed and estimated rank for two stationary targets. The results on the 
left are for one target fixed at location p x = (5, 5, 0)m. We vary the location of the other target between 
(—5, 0.01, 0)m and (—5, 30, 0)m. The range separation is 10m, so that B\f3\ = 41.47, and the g-Hankel matrix 
H has negligible entries. The results on the right are for one target at p 1 = (0.15, 5, 0)m and the location 
of the other varying between (—0.15, 0.01, 0)m to (— 0.15, 30, 0)m. The range separation is 0.3m, so that 
B \(3 1 = 1.24, and the g-Hankel matrix H is no longer negligible. We see that in spite of H being neglible or 
not, the rank of matrix C behaves essentially the same, as predicted by the asymptotic theory. This would 
be difficult to guess by just looking at the data traces displayed in Figure [12) 



19 



Fast time (seconds) xio 8 
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Figure 12: Data traces from two stationary targets. Left: the targets are at positions p 1 = (2,5,0)m 
and p 2 = (-2,15,0)m. Middle: p x = (2,5,0)m and p 2 = (-2,5,0)m. Right: p x = (0.15, 5, 0)m and 
p 2 = (-0.15,15,0)m. 



The computed and the estimated ranks grow at the same rate with the cross range of the second target, 
i.e., with |ck2 | . The growth is monotone except in the vicinity of the local minimum corresponding to the 
targets having exactly the same cross-range. In this special configuration a\ = o<2, and therefore £i = £2- 



The symbol (4.42) simplifies to 



Q{0) 



and it exceeds the threshold e||Q||oo for # in a smaller subset of (— n, 7r), than in the general case with £1 ^ £ 2 - 
The rank is defined by the size of this set, so it should have a minimum as observed in Figure [l2| 

Similar to the result in Figure [9]for a single stationary target, there is a discrepancy between the computed 
and estimated rank, due to the aperture not being large enough. This discrepancy diminishes as we increase 
n and therefore the aperture. 



4.5 Discussion 

The analysis above shows that the matrix of data traces from a moving target has much higher rank than 
that from a stationary target. The larger the speed, the higher the rank. The rank of the matrix of traces 
from a stationary target is smallest, equal to one, when the target is at the same cross-range as the reference 
point p Q . The rank increases at a linear rate with the cross-range offset from p Q . 

Comparing the results for one and two stationary targets, we see that the rank increases. The rank 
depends strongly on the cross-range offset of the targets. There is a small effect due to the separation of the 
targets in range, but it becomes negligible in the asymptotic limit n — » 00. 

Although we have not presented an analysis for more than two targets, we observe numerically that the 
rank increases as we add more and more stationary targets. The implication is that the matrix of traces 
from a stationary scene with many targets is not in general low rank. This is why the data separation with 
robust PCA should be done in successive small time windows, with each window containing the traces from 
only a few stationary targets. These traces give a matrix that is low rank, and thus can be separated from 
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Figure 13: Computed and estimated rank of matrix C for two stationary targets. One is fixed at location p x — 
(5, 5, 0)m and the location of the other varies on the line segment between (—5, 0.01, 0)m and (—5, 30, 0)m. 

the traces due to moving targets. The simulation results shown in Figure [6] illustrate this point. 

5 Numerical simulations 

We begin with the setup for the numerical simulations. Then we present three sets of results. 
5.1 Setup 



We use the GOTCHA Volumetric SAR setup described in section (4.1). The data traces are generated 



with the model (4.5). In all the simulations but the last one, the point targets are assumed identical, with 
reflectivity a q = 1. The images are obtained by computing the function (2.6) at points in the square imaging 



region of area 70 x 70 m 2 centered at p Q . The motion estimation results are obtained with the phase space 
algorithm introduced in [2]. This algorithm requires that we know the location of the target at one instant. 
We choose it at the center of the aperture, which is why there is no error in the target trajectory at s = 0. 

The principle component pursuit optimization in the robust PCA is solved with an augmented Lagrangian 
approach. It requires the computation of the top few singular values and corresponding singular vectors of 
large and sparse matrices, which we do with the software package PROPACK. 



5.2 Simulation 1 

The first simulation is for a collection of 30 stationary targets placed randomly in the imaging region, and a 
single moving target located at (0,0, 0)m at s = 0, and moving in the plane with velocity u = ^(1, l)m/s. 
The data traces are shown in Figure [5] 

The estimated moving target trajectory is shown in the left plot in Figure [l4j The blue line corresponds 
to the true trajectory. The red and green lines are the estimated trajectories with the sparse component of 
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the matrix of traces, as returned by robust PC A with and without windowing. These sparse components 
are shown in the right plots in Figure [6j The separation with robust PC A is better for the windowed traces, 
and so is the estimate of the target trajectory. This is more clear in the right plot of Figure [l4j where we 
show the error of the trajectory. 
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Figure 14: Estimation of the trajectory of a moving target in a complex scene with 30 stationary scatterers, 
placed at random in a 50 x 50 m 2 imaging region. We compare the results obtained with the sparse part of 
the data traces returned by robust PCA with and without windowing. These sparse parts are shown in the 
top right and bottom right plots in [6] 



In the left plot of Figure 15 we show the image obtained with the data traces, and with exact compensation 
of the motion of the target. The image is focused at the initial location (0,0, 0)m of the moving target, as 
expected. However, the stationary targets are out of focus, and the image appears noisy. The right plot 



of Figure 15 shows the image obtained with the sparse component of the traces, separated successfully by 
robust PCA with data windowing. The motion compensation is with the estimated velocity. We note that 
the artifacts due to the stationary targets are now removed, and the image peaks at the expected location 
(0, 0, 0)m. There are two ghost peaks, due to the error in the estimated target velocity, but they are much 
smaller than the peak at the correct location. 



5.3 Simulation 2 

The second simulation is for a scene with 20 stationary targets and two moving targets. The first moving 
target is as in the first simulation, The second one is located at (—5, 5, 0)m at s = and moves in the plane 
with velocity u = (—1, y/2)m/s. 



The data traces D r (s, t) are plotted on the left in Figure 16 The separation with robust PCA is shown in 



the middle and right plots of Figure 16 Each is normalized by the maximum of |Z) r (s, and then plotted 



on the same color scale. Note how the traces from the two moving targets are separated from the other 
traces. This simplifies the motion estimation. 
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Figure 15: Images of a scene with 30 randomly placed stationary scatterers and one moving target with 
velocity 28 m/s. Left: the image given by the original data with exact compensation of the motion of the 
target. Right: the image given by the sparse component of the data, separated from the other traces by 
robust PC A with time windowing. The images are normalized by the largest pixel value and plotted in dB. 




Figure 16: Data trace separation with robust PCA for a scene with 20 stationary targets and two moving 
targets. From left to right we show the matrix M of data traces, the low rank £, and the sparse S parts. 
In each plot we show absolute values normalized by the largest value of \D r (s,t)\. The sparse component is 
plotted in dB scale to emphasize the contrast. 



5.4 Simulation 3 



Our last simulation considers again a scene with 30 stationary targets and a single moving target. This 
target is like that in simulation one, except that its reflectivity is ten times larger than that of the stationary 
targets. 



The data separation results are in Figure 17 We show in Figure 18 the estimated target trajectory with 



the original data traces (left plot in Figure 17), and the sparse component (right plot in Figure 17). We 
obtain as before that the estimation is better after the data separation. 



The left plot in Figure 19 shows the image computed with the original SAR data traces. The image is 
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Figure 17: Data trace separation with robust PC A for a scene with 30 stationary targets and a moving one 
with reflectivity that is ten times stronger than the others. From left to right we show the matrix M of 
data traces, the low rank £, and the sparse S parts. In each plot we show absolute values normalized by the 
largest value of \D r (s,t)\. The sparse component is plotted in dB scale to emphasize the contrast. 



focused at the stationary targets, but there is a strong artifact (streak), due to the moving target. The middle 
plot in Figure 19 shows the image computed with the low rank component of the data traces, displayed in the 



middle in Figure 17 The effect of the moving target is now considerably smaller. The right plot in Figure 



19 shows the image obtained with the sparse component of the traces, with motion compensation using the 



estimated target velocity. There is no artifact due to the stationary targets and the image is focused at the 
location (0, 0, 0)m, as expected. 



6 Summary 

In this paper we consider the problem of synthetic aperture radar (SAR) imaging of complex scenes consisting 
of a few moving and many stationary targets. The SAR setup is the usual one with a single antenna mounted 
on a platform flying above the region to be imaged. With large bandwidth probing signals and long data 
acquisition trajectories, SAR can produce high resolution images of stationary scenes. However, the presence 
of moving targets may cause serious degradation of the images. When the targets have moderate speed, they 
appear displaced and blurred in the images. Fast moving targets create significant artifacts such as prominent 
streaks. 

To bring the images of the moving targets in focus we need to estimate their motion. This is necessarily 
done with small successive sub-apertures, corresponding to short acquisition times over which the target 
motion can be approximated by uniform translation. Imaging with motion estimation is difficult for at least 
the following reasons: First, the echoes from the moving targets may be overwhelmed by those from the 
stationary scenes, so the targets may be difficult to detect. Second, even if we can detect the presence of a 
moving target in a complex scene, it is difficult to estimate its motion with the existing algorithms, unless we 
use multiple receiver or transmitter antennas. The algorithms that work with the usual SAR setup assume 
that all the targets move the same way and are sensitive to the presence of strong stationary targets. Third, 
even if we detect and estimate well the target motion, when we compensate for it in the image formation 
process we may bring the stationary targets out of focus, and thus still get images with significant artifacts. 

To address these challenges, we propose a pre-processing step of the SAR data designed to separate the 
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Figure 18: Estimation of the trajectory of a moving target in a complex scene with 30 stationary scatterers, 
placed at random in a 50 x 50 m 2 imaging region. The reflectivity of the moving target is ten times stronger 
than that of the stationary ones. Left: estimated target trajectory using the original traces (green) and the 
separated traces (red). The true trajectory is in blue. Right: errors of the estimated trajectories. 
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Figure 19: Images of a scene with 30 stationary targets and a moving one with reflectivity that is ten times 
stronger than the others. Left: image obtained with the original data traces. Middle: image obtained with 
the low rank component of the data traces, returned by robust PCA. Right: image obtained with the sparse 
component of the data traces, returned by robust PCA. The motion compensation is with the estimated 
target velocity. 



stationary target echoes from those due to the moving targets. The main result of the paper is to show 
that this can be accomplished with the robust principal component analysis (PCA) algorithm complemented 
with appropriate data windowing. The robust PCA algorithm decomposes a matrix M into a low rank part 
C and a sparse part S. In our context, the matrix M is given by the pulse and range compressed echoes 
received at the SAR platform. We show with analysis and numerical simulations that the contribution of the 
stationary targets to M is a low rank matrix when we observe it in a small enough time window. Thus, we 
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may think of it as the component C of M. The contribution of a few moving targets to M is a sparse matrix 
that has higher rank, depending on the target velocity. Therefore, we expect that robust PCA separates it 
from the low rank part £, due to the stationary targets. 

We show with numerical simulations that indeed, robust PCA can accomplish such data separation. But 
the algorithm cannot be applied as a black box. It must be complemented with proper windowing of the 
pulse and range compressed SAR data in order to achieve a good separation. We present results for various 
imaging scenes containing multiple stationary targets, and one or two moving targets that may be stronger 
or weaker than the stationary ones. For weaker targets, we demonstrate that robust PCA can detect their 
faint echoes and it separates them from those due to the stationary targets. We also show motion estimation 
and imaging results with and without the data separation step, in order to demonstrate its importance in 
achieving good results. 
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Appendices 

A Single target covariance matrix 

We approximate here the function 




xcosK(t-Ar( S / ,p( S / ))))exp - — (t - Ar(s', p(s'))) 2 , 



and simplify notation as 



At" (a, s') := At(s, p(s)) ~ Ar(*', p(s')) 
At+(s, s') := At(s, p(s)) + Ar(s', p(s')) 
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We rewrite the integrand using a trigonometric identity, and completing the square 



cos[ou (t — Ar(s, p(s)))] cos[u (t — Ar(s / , exp 



B 



B 2 



-(* - Ar(s, p(,))) 2 - — (t - Ar( S , p(.))) 2 



exp 



£ 2 (Ar- (*,*')) 



cos[cjoAt (5, s')] + cos 



2uj [t + 



At+(s,s') 



x exp 



-B 2 t 



At+(s,s') 



and obtain that 



CM 



2At 



cos[cjoAt (s, s')\ exp 



B 2 (Ar-(s, S ')) 



1 

+ 2At 
2BAt 
2BAt 



exp 



5 2 (Ar-( S , S 0) 2 



4 

(it cos 



cos[o;oAr (s, s')] exp 
cos[cjo Ar~ (s, s')] exp 



B 2 (At-(s,s')Y 



B 2 {Ar-{s,s')Y 



dt exp 

exp 



„ 2 / At+(s,s') 
£ 2 It + V 



exp 



2 

-£ 2 ( * 



2£At 



£ 2 (At-( S ,s')) 



2/ , , Ar+ 



£ 2 



The last approximation is because in our regime uo Q ^> B. 



Our assumption (4.16) on the Fresnel number, and therefore on the aperture, allows us to linearize 



Ar (s, s') and obtain the Toeplitz structure stated in Proposition [T] We have by the mean value theorem 
that 



At-(s,s') = (s-s')^At(s,p(s)) 



for some s between s and s f . The derivative is given by 



j-Ar(s, — Vt(s) • (m(s) - m (s)) - u • m(s) 

as c. 



(A.l) 



(A.2) 



in terms of the unit vectors 



|r( S ")-p( S -)|' ^' |?( S -)-p | ; 
and the unit vector t(s) tangential to the flight path at r(s). We use that 



m(s) - m (s) = [/ - m (s)m (s) T ] ffij + O 



\r(s)-Po\ 



and expand the right hand side in (A.2) around s = and obtain 



d_ 
ds 



At(s,p(s)) 



2 

c„ 



f P Q (p-p Q ) 

. 'I?(0)-Po| 



u • m„ 



P (p-p ) 
' I?(0)-P„|J 



at • P u 

c L 



O 



VaR 
c n L 2 



. (A.3) 
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Therefore, 



ui At (s,s') 



2(s-s') 



P (p-p ) 

. '|?(0)-P | 



u • m n — u 



Vq(P-Po) 
\r(0)-Po\ 



with negligible error by assumption (4.16) 



8 = 



a 2 t ■ P u 
XoLV 



\xJ? <<L 



This is the result stated in Proposition [T] 



B Computation of the symbol 



The symbol is given by 



Q(0)= c i e ^> 0e(-7r,n) 



j = -oc 



with cj defined by (4.24). Thus 



r~ 00 2 



- E 



OO 



cos(jQ) cos(7j)e" 



. (iir 



45At 

j = -oo 

where we recall that 

Define Ax = £/2 and xj = jAx. Then 



2£A£ 

/ — 00 

£ [cos(n(^-7))+cos(n(^ + 7 ))]e 



E 



sin(n#) cos(7j)e 



j = -oo 

(U) 



Q{6) ~ iBAtAi ^2 



cos 



£ = B\a\As. 



2(0 -7) 



2(^ + 7) 



For G (— 7r,7r) and £ (i.e., As) small enough, we can approximate the sum with an integral 



Q(0) 



/7T 



2BAt£ 

7T 

2BAt£ 



cos I 2 ^ 7 ^ x) + cos (^—^^x 



e x dx 



(fl-7r _(!+ir 



(A.4) 



(B.l) 



C Rank estimate of large Toeplitz plus g-Hankel matrices 

We use the results from [19] to obtain the asymptotic estimate of the rank of matrix C given by equation 



(4.37) as the sum of a Toeplitz matrix T, a g-Hankel matrix H and its transpose. We need the following 
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definition: 



Definition 2 Let Q be a complex valued, measurable function defined on the interval U = (— 7r,7r). Let 
also {A n } ne ?q be a sequence of of matrices. Each matrix A n G R( n + 1 ) x ( n+1 ) ; and we denote by o-j(A n ) its 
singular values, in descending order, for j = 1, . . . , n + 1. We say that the sequence is distributed (in the 
sense of singular values) as the pair (Q, U), and write in short {A n } ~ a (Q, U), if 

n+l 

' F(\Q(0)\)d0, (C.l) 



1 _ 1 T 71 " 

.7 = 1 



/or ever?,/ F G C (1R + ). iJere C G (1R + ) is the set of continuous functions with bounded support over the 
nonnegative real numbers. 

It is shown in [19, section 4.2.2] that if {H n } ne ^ is a sequence of g-Hankel matrices H n G R( n + 1 ) x ( n + 1 ), 
then 

{H n }~ a (0,U). (C.2) 

Moreover, [19, Proposition 4.3] states that if {A n } ne ?q and {H n } ne ?q are two sequences of matrices satisfying 
{A n } ~ a (Q, U) and {H n } ~„ (0, U) then 

{A n + H n }~ a (Q,C7). (C.3) 
In our context, A n = T n are Toeplitz matrices defined by the sequence {c^jjez, and Q is the symbol defined 



by (4.26). We are interested in the distribution (in the sense of singular values) of the sequence {C n } nG ^ 
defined by 

C n = T n + H n + Hl. (C.4) 

Because the singular values of the transpose H% are the same as the singular values of H ni we obtain using 
( [021 ) and Definition ([5} that 

AW (C5) 



Thus, {(T n + # n ) + Hi) has the same distribution as {T n + i^ n } and by ( |C.3[ ), 

{T n+J ff„}~ CT (Q,[/). (C.6) 

More explicitly, 

lim 4rE%( T n + ^ + ^)) = r ^(IQWIW VF G C (K+). (C.7) 

n^oo n+l z — ^ Z7T J_ 7r 

We cannot apply directly this result to the computation of rank, because the indicator function l[ e ||Q || 00,00) 
does not have bounded support and it is not continuous. However, the result extends to such functions as 
we shown next. 
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First, let us show that the sigular values of matrices T n + H n are bounded uniformly in n. Because the 
largest singular value of a matrix is equal to its 2-norm, we obtain by the triangle inequality that 

<ri(T n + H n + Hi) = \\T n + H n + ^|| 2 < ||T n || 2 + 2\\H n \\ 2 = <ii (T n ) + 2<7i(iT n ). 

One of the results of the Szego theory for large Toeplitz matrices [T71 [3] is that the sequence of largest 
singular values {(Ji(T n )} ne ^ converges to the limit ||Q||oo- Thus, the sequence is bounded above, and we 
denote the bound by E^. For the g-Hankel matrix we can use the matrix norm inequality 

(Ti(fr„) = ||fr n || 2 < VWHJiWHjoo, 



where ||i? n ||i an d ||^n||oo are equal to the maximum of the 1-norm of the columns and rows of H ni re- 



spectively. They are defined by equation (4.39), in terms of the sequence {hj} given in equation (4.41) 



Obviously, the 1-norm of the columns and rows of H n are bounded above by the series 



3=0 



OO, 



which is convergent because hj decays exponentially with j. Therefore, 



WHJUKVh and llJ^lloo < S 



H 



and gathering the results above, we have that 

^i(T n + tf n )<X T + 2£tf. (C.8) 

Since the singular values are bounded by ( |C.8| ), in our calculation of rank we can replace the indicator 
function I^hqh^oo) ( x ) with a new function x( x ) of bounded support, satisfying 

{0 x < 5, 
(C.9) 
1 xe[5,D], 

and decaying to zero in a continuum manner for x > D. Here we simplified notation as 

5 — ellQHoo and D = E T + 2E^. 



It remains to show that result (C.7) extends to the function which has bounded support but is discontin- 
uous at x = S. 

Let us introduce the sequence of continuous functions {Xm(x)}mez+ 1 defined by 

Xm{x)= {m(x-6) + l 6-±<x<5 (CIO) 
x(x) x > S. 
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This sequence converges pointwise to x(x), as m -+ oo. We know that (C.7) holds for F — Xm- To extend 
the result to F = x> we show next that we can interchange the limits as in 



n+l 



n+1 



lim lim J^Xm (<Tj(T n + H n + H%)) = lim lim V X m (^(T n + # n + f#)) 

J = l J = l 

n+l 

= j!So^&( ^ ( T n + # n + ^ ) ) ■ 



We have that 



n+l 



n+l 



lim lim —j— Xm {*, (T n +H n + Hl))- lim — j- £ X (T„ + ff„ + H^)) 

2— ^oo n— ^oo n+l — n— )>oo n+l L — ' 
. n+l 

^ ^XT E faC^n + H n + ffj)) - X (cJ 3 (T n + tf n + !#))] 

n 7-00 n — I - 1 

J = l 
^ n+l 

i™o ^Xi E [9m + H n + Hi))] 



lim 

ra— >-oo 



lim 

ra— >-oo 



(C.ll) 



with residual g m := x™ — X that satisfies by construction 

9m{x) > 0, Vx G 
and is bounded above by the continous function 



Gm(x) < 



x < S- ±- 

m 

m(x - S) + 1 (5-^<x<(5 
-m(x - 5) + 1 S < x < S + ^ 
*>£+^. 



This gives 



< lim 

m— )-oo 



^ ^TI E M T « + + *#))] 



< lim 

m— ^oo 



n+l 



= lim + r g m (Q(0))de, 

m^oo Z7T J_ 7T 

and we can bring the limit inside the integral using the dominated convergence theorem 
lim + P Q m (Q(6)) d6=±- P lim Q m (Q(6)) d6 = 0. 

m^oo Z7T J_ 7T Z7T J _ 7T m^oo 



The last equality is because Q m — » almost everywhere. Thus, the limit (C.ll) is equal to zero, and the 
result follows. 
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